06 - Raster manipulation and analysis

What you learn today? Wrangle your rasters

In these exercises, you learn how to work with raster data in R, specifically, how to

review and recode problematic raster values
crop raster extent and mask it to a boundary specified by a polygon
reduce raster resolution
extract values from raster by points, lines or polygons
overlay multiple rasters and do raster algebra
inspect prominence at locations and sort it
sum raster values across vector shapes

Welcome to Bulgaria

In terms of data, today you will be working with an ASTER digital elevation model, an IKONOS multispectral satellite image and a prominence raster, all for the Kazanlak Valley in central Bulgaria. The images are provided by the JICA and GeoEye Foundations respectively, the prominence raster is derived from the elevation model:

Aster image contains a digital elevation model for the area.
IKONOS image is a 4-band satellite image covering ca 150 sq km of the Kazanlak Valley with red, green, blue and infrared band, captured in 2001. As it is a fairly large orthophoto, I divided the original into two halves with each half at 2 GB and placed them in ScienceData. If you want to play with the full-sized images, you can download them manually from public www.sciencedata.dk folder, or directly with with file.download() using these direct links for West and East respectively. You can also skip the large files and load a slightly reduced image of the East side of the Valley KazE.tif from your data/ folder in the course of exercise 4.
prominence is a 30m resolution raster which assigns each cell a percentage that represents the number of cells in a radius of 1500m that are lower in elevation than the given cell.

Task 1: Access raster data values

Raster data can be very big depending on the extent and resolution (grid size). In order to deal with this the raster() and brick() functions are designed to only read in the actual raster values as needed. To show that this is true, you can use the inMemory() function on an object and it will return FALSE if the values are not in memory. If you use the head() function, the raster package will read in only the values needed, not the full set of values. The raster values will be read in by default if you perform spatial analysis operations that require it or you can read in the values from a raster manually with the function getValues().

Instructions

Activate raster library
Use GDALinfo() to inspect the properties of the Aster image raster in the data folder (the file name is Aster.tif). What can you learn from this inspection about its bands, resolution, and projection?
Read in the Aster image. Review Week 02 raster loading if unsure about how.
Use the inMemory() function on the elevation object to determine if the data has been read in.
Use the head() function to look at the first few values from the elevation raster.
Use the getValues() function on the elevation object to read in all the data.
Use the hist() function to create a quick histogram of the elevation values. Note the pile of values near -9999, these should be NA (any idea why?) and we will address this later.

# library
___(raster)

# Read in the elevation layer
elevation <- ___("data/Aster.tif")

# Check if the data is in memory
___(___)

# Use head() to peak at the first few records
___(___)

# Use getValues() to read the values into a vector
vals <- ___(elevation)

# Use hist() to create a histogram of the values
___(vals)

Solution

[1] FALSE

       1     2     3     4     5     6     7     8     9    10    11    12
1  -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
2  -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
3  -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
      13    14    15    16    17    18    19    20
1  -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
2  -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
3  -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999
 [ reached getOption("max.print") -- omitted 7 rows ]

Congratulations! You now know that the raster package only reads in raster values as needed to save space in memory. You can get the values manually using the getValues() function.

Now, a new question arises, why are so many values encoded as -9999? Are we in the Mariana Trench all of sudden?

Task 2: Change values: handle missing or bad data values in rasters

There are many situations where you might need to recode raster values. You may want to change the outlier values to NA for example. In the raster package, reclassification is performed with the reclassify() function.

In the elevation raster you’ve worked with the values are meters above sea level and are supposed to range between 0 and 2500. Anything below 0 should be an NA. In this exercise you will assign any values below 0 to NA.

Instructions

Check that the package raster and the object elevation are loaded in your workspace.
Plot the elevation raster using plot().
Set up a three-column matrix with the cbind() function and values -10000, 0, NA.
Use the matrix and reclassify() to reclassify values below 0 to NA. You will need to use the argument rcl.
Plot the reclassified elevation layer to confirm there are no values below 0.

# Plot the elevation layer to see the legend values
___(elevation)

# Set up the matrix
vals <- ___

# Reclassify 
elevation_reclass <- ___(elevation, ___ = ___)

# Plot again and confirm that the legend range is 0 - 2400
___(elevation_reclass)

Solution

Good work! Knowing how to reclassify rasters will come in handy. When you get a chance you should review the help for reclassify() particularly the part that discusses how to specify the rcl argument. The three-column approach from this exercise is most common but there are other approaches.

Task 3: Crop and mask rasters on the basis of other spatial objects

Mask and crop are two similar operations that allow you to limit your raster to a specific area of interest. With crop() you limit the extent of your raster to that of your focus area. With mask() you essentially place your area of interest on top of the raster and any raster cells outside of the boundary are assigned NA values.

The Aster image covers a large area, but we are primarily interested in the areas surveyed by archaeologists, the first of which is the cluster of burial mound points and second, large survey polygons registered by local archaeologists during field survey. In this exercise you will use the extent of mounds and survey units to crop and mask the elevation raster.

On some OSs, raster package does not support sf objects. This should not happen now, but if you encounter difficulty with raster:vector interactions, it helps to know that you can convert the vector to Spatial object with, for example, as(input, "Spatial").

Instructions I - Crop raster by the extent of the mounds dataset

Load sf library
Create mounds object from the shapefile “KAZ_mounds.shp”. For a refresher on vector loading, check out Week 02 instructions.
Verify that mounds CRS matches the elevation_reclass CRS and fix with st_transform() if not.
Create a bounding box mounds_bb around the mounds with st_make_grid() following Week 03 guidelines.
Plot the mounds_bb and the mounds to see how these two objects relate.
Crop the elevation_reclass layer by the mounds object using the crop() function and create a smaller elev object
Plot the elev and mounds and the mounds_bb together.

# YOUR CODE

Solution

[1] "Modes: S4, character"                                
[2] "Attributes: < Modes: list, NULL >"                   
[3] "Attributes: < Lengths: 3, 0 >"                       
[4] "Attributes: < names for target but not for current >"
[5] "Attributes: < current is not list-like >"

Instructions II - Filter largest survey polygons

Create survey object from a shapefile called “KAZ_units.shp”
Project the survey object to match the new elev raster with st_transform()if needed.
Compute the area of the survey with st_area() and save this object as areas. What units are these?
Filter the survey units to only those above 30 ha with the filter() function. You will need to wrap areas in unclass(). Save as survey_big. Remember to have the tidyverse or dplyr library attached for filter() to work properly. Also sf might interfere so specify the dplyr::filter() if needed.

# Read in the survey object
survey <- ___(___)

# Compute the area of the survey
areas <- ___(survey)

# Filter to survey with areas > 30000
survey_big <- ___(survey, ________ > 30000)

Solution

Instructions III - Mask raster by the largest polygons

Review the plot of elev raster.
Plot the geometry of the survey_big over it.
Mask the elev layer with survey_big and save as elevation_mask. This may take a couple of seconds.
Review the plot of elevation_mask.

# Plot the elevation raster
plot(________)

# Plot the geometry of survey_big
plot(_________(survey_big))

# Mask the elev layer with survey_big and save as elevation_mask
elevation_mask <- ________(elev, mask = survey_big)

# Plot the elevation_mask -- this is a raster!
plot(elevation_mask)

Solution

# Plot the elevation raster
plot(elev)

# Plot the geometry of survey_big
plot(st_geometry(survey_big))
plot(mounds_bb, add = TRUE)

# Mask the elev layer with survey_big and save as elevation_mask
elevation_mask <- mask(elev, mask = survey_big)

# Plot the elevation_mask -- this is a raster!
plot(elevation_mask)

Nice! You ensured that layers had the same CRS, you cropped the raster, computed bounding boxes and areas for masking and filtered the data. Finally, you used mask() to mask the elevation raster to show only the large survey units.

Question:

1. What extent does the elevation raster default to after cropping by mounds vs by masking by the large polygons of survey_big?

2. Why do we not use the mounds bounding box to crop the elevation raster?

Task 4: Reduce the raster resolution (grid cell size) using aggregate()

Rasters, such as orthophotos, terrain models, or mosaiced rasters for large area, often come in resolutions far greater than you need (browse examples athttps://datafordeler.dk/dataoversigt/danmarks-hoejdemodel-dhm/overflade-praedefineret-geotiff/) . Up till now you have played with fairly small, processed rasters. Now you are getting a taste of the real thing. To reduce the computational load when running analyses, or when developing the right approach, you should use a reduced resolution raster.

The function to reduce resolution in rasters is aggregate() which, as you might guess, aggregates grid cells into larger grid cells using a user-defined function (for example, mean or max). The function used to aggregate the values is determined by the fun argument (the default being mean) and the amount of aggregation is driven by the fact argument (the default being 2).

Instructions

Load the East half of the two IKONOS images of the Kazanlak Valley from data/KazE.tif. See Welcome to Bulgaria for links to full-resolution rasters.
Plot the file you read in with the appropriate function for a multi-band raster.
Determine the raster resolution using res() and number of raster cells in the layer with ncell().
Aggregate the IKONOS image using the default for fun and with a factor of 10 and save the new raster to east_small.
Plot the new raster for comparison to the old version.
Determine the new raster resolution and the number of raster cells.

# Load the IKONOS raster
east <- ________________ 

# Plot the IKONOS raster
___(east)

# Determine the raster resolution
___(east)

# Determine the number of cells
___(east)

# Aggregate the raster
east_small <- ___(east, fact = ___)

# Plot the new east layer
___(east_small)

# Determine the new raster resolution
___(east_small)

# Determine the number of cells in the new raster
___(east_small)

Solution

# Load the IKONOS raster
east <- brick("../data/KazE.tif")

# Plot the IKONOS raster
plotRGB(east, stretch = "lin")

# Determine the raster resolution
res(east)

[1] 10 10

# Determine the number of cells
ncell(east)

[1] 2374596

# Aggregate the raster
east_small <- aggregate(east, fact = 10)

# Plot the new east layer
plotRGB(east_small, stretch = "lin")

# Determine the new raster resolution
res(east_small)

[1] 100 100

# Determine the number of cells in the new raster
ncell(east_small)

[1] 23940

Lovely job! In this example you read in a raster and then converted it to a lower resolution raster to save on the size of the object and ultimately computation power required. In this example, the raster was not too big to begin with so perhaps aggregating would not be necessary but for big rasters such as you will see in remote sensing session (next week), this can be a big help and a necessity.

Task 5: Extract raster values by location

Beyond simply masking and cropping you may want to know the actual cell values at locations of interest. You might, for example, want to know the elevation at each mound location or perhaps the mean elevation within the large survey units. This is where the extract() function comes in handy.

Usefully, and you’ll see this in a later analysis, you can feed extract() a function that will get applied to extracted cells. For example, you can use extract() to extract raster values by point, line, polygon, or neighborhood and with the fun = mean argument it will return an average cell value by neighborhood.

Instructions

Ensure mounds and elev or elevation objects are still in memory.
Use the raster function extract() to determine the elevation at each of the mounds. Assign the extracted value into a elev column in the mounds object. Beware that the extract() function exists across multiple packages, so it’s wise to use raster:: in front of it.
Look at the mounds elevation column through the plot() function as well as histogram. Do the extract() results make sense? The elevation layer values represent meters above sea level.

# Extract the elevation values at the mound locations 
________$elev  <- _______::______(elevation, mounds)

# Look at the mounds and extraction results
____(mounds[_______])
hist(____)

Solution

Great! raster::extract() is a very useful tool. It can be used with polygons as well as points and lines, and the result can be written out as a separate vector or added to an existing object as a column.

Task 6: Raster math with overlay

You will now use the elevation layer and a “prominence” layer. Prominence measures what percentage of surrounding cells are below any given location in the raster. A high percentage value thus indicates a prominent point, while a low percentage indicates a low-lying location with poor inter-visibility. Archaeologists assert that mounds have been built in elevated places with a good field of view, so let’s see whether they are mostly right :)

What you will do in this exercise is essentially: identify the most prominent locations among the registered archaeological mounds by finding areas that have both a high percentage of prominence and a high elevation. You will use two rasters and you will define an overlay function f to do the raster math with. Remember that raster algebra uses basic Boolean logic to select desirable values. Note that you can only do raster math on two rasters of the same extent, so make sure you align the elevation and prominence accordingly with the crop() function

Instructions

Make sure the elevation object exists in memory (“Aster.tif”, it is a single-band raster).
Read in the prominence raster layer from “prominence1500.tif”; it is also a single-band raster.
Plot the prominence object. Do the legend units make sense?
Specify function f, where you select elevation values greater than 400 msl and below 650m (mounds don’t appear above that elevation) and prominence values over 60%
Call overlay() on elevation and prominence. Set the fun argument to f.
If you have not cropped the two rasters to the same extent yet, you can do so inside the overlay() function

# Check in on the elevation and read in prominence layer
elevation
prominence <- ___(___)

# Plot prominence
plot(prominence)

# Function f with 2 arguments and the raster math to select specific elevation range and prominence values
f <- function(rast1, rast2) {
    rast1 _________ & rast2 ___________
}

# Align the extent of the two rasters with crop()
__________

# Do the overlay using the above define f function in the `fun` argument. 
prom_el_overlay <- ___(______, ___, fun = ___)


# Plot the result (low elevation and high prominence areas)
___(prom_el_overlay)

Solution

class      : RasterLayer 
dimensions : 590, 627, 369930  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 352483.7, 371293.7, 4712336, 4730036  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=35 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : Aster 
values     : 265, 1302  (min, max)

Congratulations! You’ve now learned to perform raster math using the raster function overlay(). You limited to areas with >400m & <650m elevation and >60% prominence, these areas should harbor the most prominent mounds (or defense-ready areas if associated with settlements) in Kazanlak.

Questions:

3. What is the actual value range in the prom_el_overlay raster?

4. What area is covered by each cell?

Task 7: Inspect the most prominent mounds

In the above exercise, we produced an elevation-prominence overlay. Mounds and other sites that sit in this overlay enjoy a strategic position vis-a-vis the rest of the valley. Which ones are they, however, and what are their real prominence values? It is hard to see at the scale of the valley and it would be good to pull the sites out.

Find the mounds that enjoy the most prominent locations as well as those that feature in the elevation-prominence overlay raster. Produce a list of the ID numbers (TRAP_Code) of all the overlay mounds and 25 most prominent mounds and plot them (expressing their prominence somehow) .

Instructions

Check that you have the mounds sf object, prominence and prom_el_overlay rasters.
Plot the prom-el-overlay and the mounds on top of each other to check visual overlap. Do the same with prominence raster and mounds.
Extract values from the elevation-prominence overlay raster and from the prominence raster at mound locations and write them to two columns:

– call the first column prom_el_overlay – name the second column prominence
Make an object of mounds that sit within these strategic high-visibility locations. How many are there?
Make an object of 25 mounds with the highest prominence values (remember arrange() and slice()?). Which TRAP_Code ids are included?
Plot both these sets of mounds using the mapview() library and compare their locations.

# YOUR CODE

Solution

Question:

5. How do the mounds with high prom_el_overlay values differ from those with high prominence?

Task 8: Practice further on Aarhus municipality

Load the Danish elevation file. Crop the elevation by Aarhus municipality boundary or by local forests. Calculate a summary of elevation among the urban shelters (shelters.rds)